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ABSTRACT 



An analysis of the effects of small impulses on Earth impacting asteroids is 
presented. The analysis is performed using a numerical routine for an exact, two-body, 
analytic solution. The solution is based on two-dimensional, two-body, Earth intersecting 
elliptical orbits. Given the asteroid eccentricity, time prior to impact and impulse 
magnitude and direction, an analysis of impulse-to-minimum-separation distance is 
generated. Impulse times prior to impact from zero to a few orbits are considered. The 
analysis is presented as three-dimensional plots of minimum separation distance as a 
function of impulse magnitude, direction, and time prior to impact. The general result is 
that for long lead times the optimal impulse occurs at the perihelia of the asteroid's orbit in 
the direction of the velocity vector, in the orbital plane. For short lead times the optimal 
impulse direction becomes more normal to the velocity vector, in the orbital plane, as the 
asteroid approaches the Earth. 
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I. 



INTRODUCTION 



The possibility that a large asteroid may impact the Earth delivering an energy in 
excess of 30 MT of TNT is a very real threat. One only has to look at the Moon in the 
night sky to gain an appreciation for the magnitude and probability of such an impact. 
How the human race deals with the threat of such an impact is of significant interest to the 
planet Earth as a whole. 

As discussed by Rather et. al. (1992), as early as 1705, when Edmund Halley 
wrote A Synopsis of the Astronomy of Comets, there has been speculation that an 
extraterrestrial object might impact the Earth, the possibility of such impacts was not 
perceived as a threat to life on Earth until the late 1940's when the Moon's craters were 
fully understood to be the result of impact events, not volcanism, and that the Earth is 
subject to the same impact hazard. Rather et. al (1992) also mention that the magnitude 
of the threat was not fully realized until the late 1970s to early 1980s. The perceived 
threat began to stand upon a solid foundation with the publication by Alvarez et. al. 

(1980) of a theory on the extinction of the dinosaurs due to the impact of a large asteroid 
or comet 65 million years ago. The theory put forth by Alverez et. al. suggests that an 
impact by a 10 km diameter asteroid at the Chicxulub site off of the Yucatan peninsula 
indirectly caused the extinction of 60% of all life on the Earth, including the dinosaurs. 

Since the awareness of the possibility of an asteroid impacting the Earth began 
developing in the mid 20 th century, there have been several "near misses" recorded. 
Perhaps the most spectacular near miss was a large fireball created by an object racing 
across the daytime sky in a northerly direction that entered the Earth's atmosphere above 
Utah in the United States and exited above Alberta, Canada on the 10 th of August, 1972. 
The object observed was determined to be an asteroid upwards of 30 meters in diameter 
as reported in Sky and Telescope Magazine (1972). Had this asteroid's trajectory been 
ever so slightly different, mankind would have had its first opportunity to observe a large 
object impact the Earth. Such an impacting object would carve out a crater 200 to 500 
meters in diameter. Since then, several other asteroids and comets have been detected 
passing by the Earth at distances less than a few hundred thousand kilometers. Apollo 
Asteroid 1989FC, referenced by the AIAA Space Systems Technical Committee; Asteroid 
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1991 BA, reference by Scotti et. al. (1991); and Asteroid 1996JA1 referenced by Jaroff 
(1996) are three asteroids discovered recently passing very close to Earth. On an 
astronomical scale, these close approaches are essentially impacts. 

In April of 1990 the American Institute of Aeronautics and Astronautics issued the 
position paper entitled Dealing with the Threat of an Asteroid Striking the Earth that 
briefly described the implications of an asteroid impact. The AIAA found that "Earth- 
orbit-crossing asteroids clearly present a danger to the Earth and its inhabitants." It was 
recommended that a "systematic and open program" for detection of Earth crossing 
asteroids be established as well as a study to "define systems which can deflect or destroy, 
or significantly alter the orbits of, asteroids predicted to impact the Earth." A few search 
programs existed prior to the position paper and a few others have begun subsequently. 

As a result of the awareness of the possibility of an asteroid or comet impacting 
the Earth, the National Aeronautics and Space Administration has sponsored two 
workshops to study the fundamentals of the impact and impact mitigation problem. The 
first workshop is summarized in the report entitled "The Spaceguard Survey: Report of 
the NASA-International Near-Earth-Object Detection Workshop." The second workshop 
is summarized in the report entitled "Near-Earth-Object Interception Workshop." The 
concentration of the workshops have been related to assessing the magnitude of the threat, 
impact effects and hazards to the Earth, as well as the political implications of developing 
an impact mitigation capability. Several books have been published on the matter, both 
technical and non-technical, and more than one Hollywood movie has been based on the 
subject. Two spacecraft exploration missions. Near Earth Asteroid Rendevouz (NEAR) 
and Clementine, have included intercepts of asteroids as a major part of their mission to 
study the nature of asteroids. However, little "astrodynamical analysis" has been 
performed on the feasibility of the impact mitigation problem. 

The astrodynamical analysis of feasibility is where this thesis is intended to fit into 
the larger problem. Presented is an analysis of the impact and mitigation problem based 
on a two-dimensional and two-body analysis. It is intended to be a first order 
approximation for the solution of a larger problem. However, it is a more rigorous 
treatment of the astrodynamics of the hazard than previous published analyses, such as 
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Ahrens and Harris (1994). Included in the following analysis is the periodic nature of the 
problem as well as the near term effects. The analysis is not concerned with object 
detection or orbit prediction, but instead centers on how impulses applied to an asteroid at 
various points on the asteroid's orbit affect the outcome when there is a presumption of 
collision otherwise. Mission design for mitigation is not a goal of the following analysis; 
however, the analysis tool presented may be utilized in determining a first order estimate 
for optimizing the time and position of asteroid intercept for impact mitigation. 

Presented first is an astronomical development of the existence of an asteroid 
impact problem from the origins of the solar system to the record of past impacts on the 
Earth followed by a brief description of two-body orbits. The problem is then presented 
with governing assumptions and a method of solution. The solution method is then 
assessed and validated followed by an analysis of an asteroid impact scenario with a 
discussion of the results. The analysis method is then applied to the asteroid Toutatis 
which will make several close approaches to the Earth over the next decade. 
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II. SOLAR SYSTEM OBJECT ORIGINS 



A. SOLAR SYSTEM FORMATION 

Current theories of formation of the Solar System stem from the post initial 
expansion universe environment of gas, dust, radiation, and magnetic fields in a 
nonuniform distribution. Bouyed by the interplay of gravitational, magnetic, and pressure 
forces local mass concentrations began to form. This interplay of forces on the non- 
uniform environment sets up initial angular momentum conditions for very large three- 
dimensional structures. These large rotating structures began to coalesce into accretion 
disks around central masses. These central masses eventually reached critical mass for 
fusion to occur and stars emerged. We call our local star the Sun. 

B. PLANETARY FORMATION 

In a similar fashion to the stellar formation, the accretion disk around the Sun 
provided the environment for smaller mass concentrations and accretion disks to form thus 
producing small scale structures called planetismals composed of the basic chemical 
elements in varying quantities. Gravitational attraction and relative motion of these 
planetismals caused them to collide with one another and cohesive forces enabled some to 
remain attached forming larger structures. After enough of these interactions occurred, 
the planets began to form. Unlike for the stellar conditions, the planets do not possess the 
critical mass to initiate and sustain a fusion reaction. This allows for large scale assembly 
of solid, aqueous and gaseous structures to take place in quantities and composition 
proportional the relative percentages of chemical elements present in the planetismals. 
These structures are more familiar on the Earth as the crust, the oceans and the 
atmosphere. Analysis of the elements present from the impact delivery mechanism and 
current known compositions of asteroids and comets suggests that this was the mechanism 
of organic and non-organic material delivery that formed the Earth as proposed by Chyba 
et. al. (1994). 

C. ASTEROID AND COMET FORMATION 

However, not all of the accretion disk surrounding the Sun coalesced into either 
the Sun, the Planets or their satellites. The remainder of the planetismals continue to be 
dispersed throughout the Solar System in the form of asteroids and comets. The asteroids 
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being located primarily within the inner Solar System, and the comets existing in the outer 
Solar System and beyond. 

The asteroids are mainly concentrated in the asteroid belt located between the 
orbits of Mars and Jupiter. The remainder of the asteroids are dispersed throughout the 
solar system in elliptical orbits lying mainly inside the orbit of Jupiter. There are several 
theories of how the asteroid belt came to exist. These ideas include the destruction of a 
planet, or the inability of a planet to form, due to the tidal forces of Jupiter. There are 
several theories accounting for the formation of the asteroids not located in the asteroid 
belt. These ideas range from planetismals that have never collided and attached 
themselves to another planetary body, to cast off remnants of massive collisions of 
planetary bodies with very large planetismals, to objects from the asteroid belt that may 
have been perturbed by a passing object into a smaller orbit. All of the ideas have some 
amount of merit that give them validity. 

Comets are believed to originate from the Oort Cloud of cometary material 
orbiting the Sim at a distance of some 50,000AU. From this cloud, comets are believed to 
be injected into the solar system by orbital perturbations due to the gravitational field of 
passing stars. Once injected into the solar system, gravitational encounters with the Sun 
and Jupiter may further perturb the comets’ orbit and either “capture” the comet, so it 
remains within the solar system, or “assist” the comet on its way to interstellar space. 
Additionally, there is believed to exist a band of icy objects that extends from the orbit of 
Neptune at 30 AU out to as much as 100 AU. These objects are said to be located in the 
Kuiper Belt, so named after Gerard Kuiper, who first proposed their existence in 1 95 1 . 
The objects that form the Kuiper belt range in size and orbital characteristics to the extent 
that it is believed that Pluto may actually be a member of this group as discussed by Jewitt 
and Luu (1995). 

D. ASTEROIDS 

1. Location 

The asteroids are of particular interest as potential impact hazards in that they have 
a greater mass density than comets and are more likely to reach the surface of the Earth in 
a given impact scenario. The majority of the asteroids are located in the asteroid belt 
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between the orbits of Mars and Jupiter. Very few asteroids have been detected beyond the 
orbit of Jupiter. This lack of detection may only be an effect of the limited capability of 
current detection sensors. However, a significant number of asteroids in orbits smaller 
than that of the typical asteroid belt object have been identified. These are the objects of 
primary concern for the problem of mitigation. 

2. Quantities 

Literally thousands of asteroids have been observed and identified orbiting the Sun. 
Of those, more than 300 are considered near Earth asteroids (NEA’s) and pose a threat as 
a potential impacting object. Of greater concern are the subset of NEA’s dubbed Earth 
crossing asteroids (ECA’s) which are currently about 200 in number. The orbits of the 
ECA’s are such that they allow for the possibility of impact with the Earth at some future 
date. These ECA’s range in size from 10 km down to 0. 1 km, which is currently the limit 
of detection of asteroids by Earth based sensors. Estimates for ECA’s near and above 10 
km in diameter indicate that all of the asteroids have been identified and that for the ECA’s 
of 1 km in diameter or less, the identified asteroids represent only about 10% of those that 
are believed to exist, as stated by Grieve and Shoemaker (1994). 

3. Classification 

Classification of the ECA's are determined with respect to the Earth's orbital 
extrema. The Earth's perihelion and aphelion distances are 0.9833 and 1.0167 AU 
respectively. The orbits of the ECA’s all have perihelia less than the aphelion of the Earth 
orbit and aphelia greater than the perihelion of the Earth orbit. The ECA’s have been 
subdivided into three classes, the Atens, Apollos and Amors, based on their orbital 
characteristics as discussed by Rabinowitz et. al. (1994). Table 1 summarizes the 
distinction between the three classes. Figures 1 , 2 and 3 depict typical orbits of each of 
the three classes. 
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Class Name 


Semi-Major Axis 
a (AU) 


Perihelion Distance 

q (AU) 


Aphelion Distance 

Q (AU) 


Aten 


< 1 


- 


> 0.9833 


Apollo 


> 1 


<1.0167 


- 


Amor 


- 


1.0167 <q< 1.3 


- 


Tal 


t)le 1 . ECA Classes 



2062 Aten, Semi-major axis = 0.9666 AU, Eccentricity = 0.1826 




Figure 1. Typical Aten Orbit 
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Y (All) 



4179 Toutatis, Semi-major axis = 2.5154 AU, Eccentricity = 0.6361 




Figure 2. Typical Apollo Orbit 



433 Eros, Semi-major axis = 1.4582 AU, Eccentricity = 0.2229 




Figure 3. Typical Amor Orbit 
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4. Physical Properties 

The physical properties of the asteroids are generally that of “rocky”, irregularly 
shaped spinning objects. Estimates of asteroid densities range from a more cometary 
density of 2xl0 3 kg m' 3 to a dense asteroid of 5xl0 3 kg m' 3 with a mean density of about 
3x1 0 3 kg m 3 . Asteroids have been observed that are composed of a single solid mass as 
well as multiple mass centers either physically connected or gravitationally bound at a 
contact surface. As stated by Winters (1996), some of the asteroids may actually be 
aggregates of numerous smaller bodies that are gravitationally bound together. This 
hypotheses is further supported by analysis of object spin motion. It appears the many 
asteroids spin at angular rates sufficiently slow as to permit gravitational binding. The 
case of the comet Shoemaker-Levy 9 supports these theories in that it was gravitationally 
separated by tidal forces from a previous passage of Jupiter prior to its 1994 impact. 

E. COMETS 

1. Location 

Early in the 20 th century, Jan Oort hypothesized that virtually all of the comets 
originate in a cloud of cometary matter beginning some 50,000AU from the Sun. The so 
called Oort cloud is assumed to be a uniformly distributed spherical shell that surrounds 
the solar system and extends to approximately half of the distance to the nearest star, 
Alpha Centauri. From this cloud, comets are injected into the solar system by orbital 
perturbations of passing stars. Once inside the solar system the comets may be further 
perturbed by the planets. The planetary perturbations may “capture” the comet in the 
interior of the solar system or it may assist the comet on its way ejecting it from the solar 
system for all time. 

2. Quantities 

The number of comets in the Oort cloud is believed to be diminis hing ^ however the 
remaining number is believed to be on the order of 10 12 . The quantity of comets that exist 
within the solar system is far smaller, only 1 00 or fewer have been identified. 

3. Classification 

The “captured” comets are classified as periodic comets and are subdivided into 
two groups as summarized by Shoemaker et. al. (1994). Members of “Jupiter’s Family” 
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have aphelia close to Jupiter’s mean orbital distance from the Sun and have periods of less 
than 20 years. Members of the “Halley family” are the so called long period comets and 
have orbital periods greater than 20 years but less than 200 years. 

4. Physical Properties 

Comets are primarily composed of “rock” and “ice.” As such they have been 
sometimes called “dirty snowballs” or “icy dirtballs” depending on their relative 
composition. The icy material is believed to be composed of water, methane or ammonia. 
The rocky material is a variety of carbonaceous substances. Their mean density is less 
than that of the asteroids and is estimated to be about 1000 to 2000 kg in 3 . Asa comet 
“bums” off its icy material from repeated encounters with the Sun, the nucleus may in fact 
become an asteroid of a somewhat smaller dimension. 
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III. IMPACT RECORD 



Occasionally an asteroid or a comet’s orbit is such that it actually hits another 
object within solar system such as the event, widely celebrated in the media, of comet 
Shoemaker-Levy 9 impacting the planet Jupiter on the 16 th of July, 1994. The very same 
mechanism that caused the solar system, in particular the Earth, to form and sustain life 
now threatens the very same life. 

A. LUNAR IMPACT RECORD 

One theory of the origin of the Moon, as mentioned by Chapman and Morrison 
(1994) describes it forming as the result of a Mars size object impacting the Earth early in 
its existence . The ejecta from the event is believed to have behaved in such a fashion as 
to remain in orbit and coalesce into what is now the Moon. Further evidence of such 
impacts is the cratered face of the lunar surface. The Moon displays literally millions of 
impact sites. Since the Moon has no atmosphere or large scale geologic processes there is 
no erosion of the impact record, unless by subsequent impact the previous impact 
structure is destroyed. Thus, the Moon serves as a reminder for all time of the nature of 
the impact hazard. 

B. TERRESTRIAL IMPACT RECORD 

On our planet Earth, the atmosphere, oceans, volcanism and plate tectonics tend to 
erode past impact sites. As discussed by Grieve and Shoemaker (1994) there are currently 
about 140 known impact craters around the world. The locations of these sites are 
displayed in Figure 4. 
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Figure 4. Known Impact Site Locations from Grieve and Shoemaker (1994) 



These impact craters range in size from 1.5 km in up to 200 km in diameter. 
Perhaps the best example of a classic impact crater is the 1 .5 km diameter Barringer Crater 
located in Arizona, see Figure 5. Barringer Crater is believed to be the most recent impact 
site on Earth having formed around 50,000 years ago. 




Figure 5. Barringer Crater from Grieve and Shoemaker (1994) 
However, the most recent impact event to have resulted in surface damage is 
believed to be the Tunguska event which took place over northern Siberia in 1908. The 
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Tunguska event is believed to be the result of a 60 m diameter comet or light asteroid that 
exploded 10 km up in the atmosphere releasing some 30 MT of energy leveling 2500 
square kilometers of forest. The Tunguska event is characteristic of a small scale impact 
occurrance. On the large scale end of the impact spectrum lies the K/T impact event so 
called by its occurrence in time at the boundary of the Cretaceous and Tertiary periods. It 
is believed that an enormous asteroid some 10 kilometers in diameter created the 200 km 
diameter Chicxulub impact site beneath the Gulf of Mexico off the coast of the Yucutan 
peninsula and is responsible for the extinction of 60% of the living species on Earth 65 
million years ago. 

There have been numerous recent close calls of asteroids impacting the Earth. 
Jaroff (1996) describes a near impact as recently as June of 1996 where an object roughly 
600m in diameter passed within 450,000 kilometers (about 70 Earth radii) of the Earth. 

C. OTHER PLANETARY IMPACTS 

Numerous impact sites have been observed by planetary exploration spacecraft that 
have been sent to Venus and Mars. As mentioned above, Comet Shoemaker-Levy 9 
impacted the planet Jupiter in July of 1994. 

D. IMPACT SCALE 

The energies released by the impact of asteroids on the Earth are quite large. It is 
conventional to express the impact energy in terms of megatons of TNT (MT), where 
1MT = 4.2x10 15 J. With the assumption of a mean asteroid density of 3000 kg m' 3 a size 
versus impact velocity relation may be made for a given impact energy. Figure 6 displays 
an estimate of asteroid mass versus impact velocity for a range of impact energies from 
tens of megatons of TNT up to a petaton (10 15 ) of TNT. The diameter estimate assumes 
an effective spherical radius corresponding to the mass for the same impact energy and 
velocity. 
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Figure 6. Impact Scale 

Estimates of impact energy and terrestrial devastation have lead to the 
classification of impacts as local, regional, and global events as described in Morrison et. 
al. (1994). The distinction between these is somewhat blurred, however it is proposed 
that local events correspond to impact energies in the vicinity of 30 MT or less which will 
affect approximately 0.001% of the Earth’s surface area or about the size of a large 
metropolitan area. Regional impact events are those that release energy in the vicinity of 
300 MT to 3x1 0 4 MT and affect about 0.1% of the Earth’s surface or about the size of a 
large state. Global events are considered to be impacts that release energies near and 
above 3xl0 s MT affecting approximately 10% of the Earth’s surface area or about the size 
of a large country. These large events may cause such disruption of the ecological 
environment by particulate injection into the atmosphere that greater that 25% of the 
Earth’s population may be eliminated. Figure 6 also shows several representative impact 
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scenarios for a typical impact velocity of 20 km s' 1 . The Tunguska event is depicted by the 
symbol. Tunguska represents a nearly classical local impact event of an object some 60 
m in diameter. The recent May 1 996 near miss is depicted by the '+' symbol and 
represents a regional event for an object some 600 m in diameter. A global impact is 
indicated by the 'o’ symbol and corresponds to an asteroid some 1 .5 km in diameter. The 
'©' symbol represents the K/T event that created the Chicxulub impact site and 
corresponds to a 10 km diameter object. It is noted that the K/T event is well above the 
threshold for global catastrophe. 
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IV. ORBITS 



A. CONIC SECTIONS 

In a simple inverse square gravitational field, orbits take the shape of conic 
sections. For an object that lacks sufficient energy to escape the gravitational attraction of 
the central body, the resultant orbit will take the shape of an ellipse. For an object with 
sufficient energy to escape the gravitational attraction of the central body, the resulting 
orbit will take the shape of an hyperbola. The transition from an elliptical to hyperbolic 
orbit occurs along an “escape” trajectory shaped as a parabola. Straight line orbits that 
lead to collision of the orbiting body with its mass center are special cases of elliptical, 
parabolic, and hyperbolic orbits. 



B. COORDINATE FRAMES 



1. Three-dimensional 

To define the physical problem in time and space a reference frame needs to be 
established. For the general three-dimensional case a Cartesian Sun-centered inertial, or 
Heliocentric, coordinate frame is chosen. Figure 7 displays the orientation of the 
Heliocentric coordinate frame. 




Figure 7. Heliocentric Coordinate System from Bate et. al. (1971) 
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The direction in the ecliptic plane from the Sun to the First point of Aries serves as 
the primary coordinate direction in the Heliocentric frame. The ecliptic normal serves as 
the “vertical” coordinate and is defined as positive in the northern half of the celestial 
sphere. The last coordinate direction is defined by taking the cross product of the 
previous two coordinate directions. 

2. Two-dimensional 

To simplify matters for the present analysis, a two-dimensional planar perifocal 
coordinate frame is chosen. The principal axis defining a two-dimensional elliptical orbit is 
the direction toward periapsis from the primary focus. It is from this primary axis that the 
true anomaly is measured in a counterclockwise sense. The secondary axis is normal to 
the primary axis in a right handed sense. The resultant coordinate system is shown in 
Figure 8. 



Q 




Figure 8. Perifocal Coordinate System 
C. ORBITAL ELEMENTS 

1. Three-dimensional 

In the general three-dimensional case six orbital elements are required to define a 
particular location in space of an object in orbit. Those orbital elements are the semi- 
major axis (a), eccentricity (e), inclination to the ecliptic (i), longitude of ascending node 
(Q), argument of periapsis (©) and the time of periapsis passage (T). See Figure 9. 
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Figure 9. Three-dimensional Coordinates from Bate et. al. (1971) 

2. Two-dimensional 

In the perifocal coordinate system, only the semimajor axis (a), eccentricity (e), 
and true anomaly (v) are required to fix a position on an orbit. (Note the argument of 
periapsis is defined to be zero in this coordinate frame.). Other parameters of concern are 
the perihelion distance (q) and the aphelion distance (Q). Figure 10 displays the two- 
dimensional perifocal coordinate system. 




D. KEPLER’S SECOND LAW 

Motion of an object about the primary focus in an elliptical orbit is governed by 
Kepler’s Second Law: the line joining the planet to the Sun sweeps out equal areas in 
equal times. This relation increases the difficulty in the problem solution in that there 
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exists a trancendental relationship between time and position. This relationship is called 
Kepler’s Equation and takes the form: 

E -esinE 



_ e + cos v 

where cos E = . 

1 + ecosv 



( 2 ) 
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V. 



PROBLEM FORMULATION AND SOLUTION 



A. STATEMENT 

Given an impending asteroid impact with the Earth, one would like to know if 
there is an optimum point on the asteroid’s orbit where an impulse may be applied to yield 
the greatest achievable separation distance at the closest point of approach for a fixed 
impulse. Additionally, it is desirable to determine if is there an optimum direction 
associated with the applied impulse that further increases the separation distance. 

B. ASSUMPTIONS 

In order to proceed with an approximate solution method, simplifying assumptions 
were made. The first assumption was that of two-body motion. That is, the Sun is the 
mass center about which the asteroid and the Earth orbit. Furthermore, the asteroid and 
Earth do not gravitationally interact with each other. It was also assumed that there were 
no external perturbing effects on the orbits due to non-gravitational forces other than the 
applied perturbing impulse. All orbits were assumed to be coplanar, which yields a two- 
dimensional problem. The perturbing impulse is assumed to occur instantaneously. The 
asteroid is assumed to be one of the near Earth objects and hence in an elliptical orbit 
around the Sun. Hyperbolic and parabolic orbits were not considered for this analysis. 
Finally, it was assumed that the Earth is in a perfectly circular orbit at 1 AU. 

C. TEMPORAL CONSIDERATION 

In determining the separation of a NEO from the Earth, time becomes the 
dominant factor in solving the problem. The relative phase of each of the orbiting bodies 
determines whether the bodies will collide. Hence, the Earth-to-NEO separation distance 
is the quantity of concern. Changing the orbital elements becomes secondary to changing 
the asteroid's orbital phase with respect to the Earth. 

D. METHOD OF SOLUTION 

Solution to the above problem may be achieved by use of a numerical simulation 
scheme where the orbital equations of motion are integrated from some initial condition 
forward in time. However, such a simulation is time consuming (on the order of several 
minutes per impact scenario) and therefore limits the scope of any analysis. The 
assumptions stated above allow for use of analytic elliptical orbit equations. Building a 
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solution based on these analytic two-body equations offers a greatly reduced time for 
solution (on the order of a couple seconds per impact scenario) and therefore a greater 
scope of analysis may be performed. A Mathworks MATLAB model was constructed to 
numerically execute the following method. A numerical integration simulation was also 
developed in order to validate the analytic method. The solution description below is 
quite general. For a detailed solution description see Appendix A. The MATLAB script 
that corresponds to the solution method is displayed in Appendix B. 

1. Geometry 

A general description of an elliptical orbit intersecting a circular orbit is used in the 
problem solution. This description suffices for all planar intercept scenarios. A rotation of 
the perifocal coordinates may be required to bring the model in alignment with the 
particular problem, but the relative geometry remains unchanged. Figure 1 1 demonstrates 
this equivalence. 




Figure 1 1 . Equivalent Impact Scenarios 

To uniquely fix an intercept scenario with the above assumptions, only the impact 
true anomaly and orbital eccentricity need to be defined. In the case of planar , orbits, the 
perihelion direction, as defined by original asteroid elliptical orbit, defines the principal axis 
from which the impact true anomaly is measured. Implicit in the above impact location 
description is the assumption of an Earth orbital radius of 1 AU. 
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2. Solution Flowchart 

Figure 12 depicts the steps in the problem solution method. 




Figure 12. Solution Flowchart 
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3. Given Conditions 

For a given impact scenario, it is assumed the impact true anomaly, asteroid orbital 
eccentricity are known and the desired impulse time prior to impact is specified. It should 
be noted here that specifying impact true anomaly and orbital eccentricity fixes the semi- 
major axis and hence, the orbital period. If either parameter is changed, the semi-major 
axis changes. That is, as impact true anomaly increases from 0 to n the semi-major axis 
decreases, thus decreasing the orbital period. 



4. Unperturbed Orbital Elements 

From the given conditions and the known Earth orbit the perihelion distance, semi- 
major axis, and orbital period for the unperturbed orbit may be determined by: 

rJl + ecos v.) 

r = — — 

p 1 + e 




and P = — , 
n 

where n = 

5. Impact Condition 

The time from perihelion of the impact with respect to the unperturbed orbit may 
be determined by Kepler’s Equation as given in Equations (1) and (2). 

6. Initial Impulse Condition 

The time from perihelion of the impulse with respect to the unperturbed orbit may 
be determined by Figure 13 and the equation: 

(Impulse Timejpenheiion = (Impact Time)p er jheiion * (Impulse Timejimpact Time- 
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Figure 13. Impulse to Impact Time Relation 



The true anomaly at impulse may be determined by inverting Kepler’s Equation 
such that true anomaly becomes a function of time. This approach requires an iteration on 
eccentric anomaly. 

The impulse position, r , and velocity, v , may now be found from: 



7. Perturbation 

An orbital perturbation, Av , is applied to the asteroid at the impulse position r . 
This yields a perturbed orbital velocity v + Av . 

8. Perturbed Orbital Elements 

From the impulse position r and perturbed velocity v + Av the perturbed orbital 
elements may be determined from: 




(3) 



1 + ecos v ’ 

r = (r cos v)x + (r sin v)y , 



(4) 




where p = a(l - e 2 j. 



h = r x (v + Av) , 

h 2 




e = — 




and Equation (4). 
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9. Perturbed Impulse Condition 

Substituting the perturbed orbital elements into Equation (3) enables the true 
anomaly at impulse to be determined with respect to the perturbed orbit. The time from 
perihelion of impulse with respect to the perturbed orbit may be found from Kepler’s 
Equation. 

10. Orbital Positions at Impact Time 

To determine the separation distance of the perturbed asteroid from the Earth, the 
positions, and times, of the unperturbed asteroid orbit about the impact point must be 
mapped onto the corresponding positions and times of the perturbed asteroid orbit. That 
is, a one-for-one correlation between where the asteroid would have been if not for the 
perturbation and where the asteroid is due to the perturbation must be developed. This is 
the key step in determining the effect of the impulse perturbation. 

The true anomalies and time from perihelion of the perturbed and unperturbed 
orbits are related by Kepler’s Equation and not a simple function. The approach used to 
achieve the required mapping is as follows: 

a) Unperturbed Orbital Positions 

Develop an evenly spaced window, or mesh, of time about the impact 
position wide enough to include the perturbed orbits’ minimum Earth separation. (A first 
estimate of the required width of this mesh is achieved by a numerical simulation. From 
the numerical simulation and model development an interval width of ± 1.5 x Av x 
Impulse TirneimpactTime was used to provide a sufficiently wide mesh to obtain a solution 
without excessive computation time.) By inverting Kepler’s Equation, the true anomalies 
of the mesh points may be determined. From the true anomalies of the mesh points the 
orbital positions may be found from Equations (3) and (4). 

b) Perturbed Orbital Positions 

From the relationship of the time of impact (known) with the time of 
impulse (also known) the center of the mesh may be determined for the perturbed orbit. 
Again, the true anomalies of the mesh points and the orbital positions for the perturbed 
orbit may be determined by use of Kepler’s Equation and Equations (3) and (4). This now 
yields the asteroid orbital positions due to the perturbation. 
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c) Earth Orbital Positions 

In the same fashion, the Earth orbital positions corresponding to the mesh 
points may be determined. However, since a circular Earth orbit was chosen, it is easier to 
relate the mesh times to orbital positions by use of the Earth’s orbital mean motion. 

11. Earth to Asteroid Separation Distance 

With the perturbed asteroid’s orbital positions and the corresponding Earth orbital 
positions known, it is a simple matter to determine the Earth to asteroid separation 
distance at each of the mesh points from: 

R i/sep = ^( x i / r- x i /E ) 2 +(yi.r-yi,E ) 2 • 

12. Minimum Separation 

From the above orbital separation distances, the minimum may be determined. It is 
useful to express this separation in terms of Earth radii. This enables a quick evaluation of 
whether a sufficient separation was achieved to cause a “miss”. As with the analysis 
performed by Ahrens and Harris (1994) the resultant separation distances scale linearly 
with the applied impulse magnitude. 

13. Sample Model Output 

A sample of the MATLAB model output is shown in Figures 14 and 15. The 
impact scenario is for the case of an impact occurring at a true anomaly of 30°, an asteroid 
orbital eccentricity of !4, and a perturbing impulse of 1 m s' 1 in the direction of the velocity 
vector occurring 0.47 asteroid orbits prior to impact. Figure 14 displays the initial 
conditions of the scenario at the time of impulse. The impact position and the asteroid and 
Earth positions at the time of impulse are also displayed. Times prior to impact are 
displayed on the asteroid orbit in tenths of an orbital period. 
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Impact True Anomaly = 30 Orbital Eccentricity = 0.5 




Figure 14. Initial Conditions 



Figure 1 5 shows the effect of the impulse on the asteroid near the impact point. 
The unperturbed asteroid position is shown achieving impact conditions at zero Earth 
radii. The perturbed asteroid trajectory is shown approaching an impact condition, but 
instead reaches a minimum separation (indicated by the 'x') and then increases in 
separation. 



Perturbed v. Unperturbed Asteroid to Earth Separation 




Figure 15. Perturbed Asteroid to Earth Separation 
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The minimum separation indicated in Figure 15 is the final output of the model. 
The figures were generated from the model working variables. This enables the model to 
be incorporated into a controlling routine allowing for repeated simulations that sweep 
over ranges of the input conditions. 
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VI. PROGRAM VALIDATION 



Prior to making an analysis of various impact scenarios, the solution method and 
MATLAB model had to be validated. This validation was achieved by use of a simple 
numerical integration simulation and by comparison to approximate analytic solutions of 
nearly circular orbits. 

A. NUMERICAL SIMULATION 

The numerical integration simulation was developed usingthe Mathworks 
MATLAB and SIMULINK numerical processors. The SIMULINK diagram and 
MATLAB script files that support the SIMULINK model are displayed in Appendix C. 
The two-body equation of motion integrated by the model is: 



r = 




Expressed as a two-dimensional system of linear differential equations in perifocal 
coordinate form: 




x = x. 




y = y. 



Numerous orbital scenarios were run to ensure the numerical integration was 
performing correctly. A fourth order Runge-Kutta integration scheme was used for the 
simulation acting on a second order differential equation. This yields a solution that is 
analytically exact and accurate to the numerical precision of the computer microprocessor. 
For the present analysis, the computer microprocessor was a 100 MHz, 32 bit, Intel 
Pentium. For the cases chosen for orbital modeling verification, integrating around one 
orbit resulted in an ending position the same as the starting position with a relative error of 
2x1 O' 16 (3x1 O' 7 km error at 1.5xl0 8 km, 1 AU). 
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Once the numerical integration orbit models were validated, numerous perturbation 
simulations were made to provide a reference data base for the analytic method. It was 
found that the analytic method numerical model was in excellent agreement with the 
numerical simulations. Differences arose only from the difference in mesh size between 
the two solution methods, with the analytic solution having a finer mesh interval. 

B. APPROXIMATE ANALYTIC SOLUTIONS 

For circular orbits, Ahrens and Harris (1994) performed an approximate analysis 
that yields expressions for maximum orbital separation due to impulses applied both 
normal to and parallel to the orbital velocity vector. For the case where the impulse is 
applied normal to velocity vector, the maximum separation is found half an orbital period 
later with a magnitude 5^ = 2AvP / n . This case is shown in Figure 16. 



|« § max * 2AvP/n 




Figure 16. Impulse Normal to Velocity Vector 



For the case where the impulse is applied parallel to velocity vector the separation 
a full orbital period later is found to be 5 = 3AvP per orbit. That is, for this case the 
separation increases by 8 for every orbit after the single impulse. This case is shown in 
Figure 17. 
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Av 



8 « 3AvP 

per orbit 



Figure 17. Impulse Parallel to Velocity Vector 

Evaluation of these scenarios by numerical simulation verifies the above 
approximate solution. The separations between the perturbed and unperturbed orbits at 
the described locations does indeed agree very well with the approximate analytic solution. 
Further, evaluation of the above scenarios by the analytic method numerical model using 
an asteroid orbital eccentricity of 10' 5 (nearly circular) and an impulse of 1 m s' 1 yields 
virtually the same result as the approximate analytic solution. Table 2 summarizes the 
performance of the three solution methods for the circular orbit case. 



Impulse Direction 


Approximate 
Analytic Solution 


Analytic Numerical 
Method 


Numerical 

Simulation 


Normal to v 


3AvP = 14.8 


14.8 


14.8 


Parallel to v 


2AvP/ rc=3.15 


3.14 


3.15 



(Separation in Earth Radii) 

Table 2. Validation of Solution Method 
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VII. ANALYSIS AND RESULTS 

A. ANALYSIS 

The preceding analytic numerical method for solution of the Earth to asteroid 
separation problem was incorporated as a function into a routine that sweeps over impulse 
direction and time of impulse. In this manner hundreds of solutions may be generated in a 
few minutes. Appendix D contains numerous analyses for impact scenarios that may 
occur around the Earth’s orbit. 

1. Long Time Response Behavior 

For a perturbing impulse time well prior to impact the observed behavior resulting 
from the orbital dynamics is not linear. The Earth-to-asteroid separation achieved by an 
impulse is strongly dependent on the location of the impulse on the orbit as well as the 
direction of the impulse with respect to the orbital velocity. The typical behavior of an 
impulse is depicted in Figure 18. Each point on this plot represents the minimum 
separation point depicted on Figure 15. In Figure 1 8 the vertical axis represents the 
minimum Earth-to-asteroid separation (in Earth radii) that results from an impulse of 1 m 
s' 1 . The two horizontal axes represent the direction of the impulse and the time prior to 
impact of the impulse. The axis representing direction of impulse is measured from 0° to 
360° with respect to the forward direction of the velocity vector. The sense of the 
direction is that impulse directions between 0° and 1 80° point inward to the orbit and 
directions from 180° to 360° point outward from the orbit. The remaining axis represents 
impulse time prior to impact as a fraction of the asteroid's original orbital period (e.g., t/P 
= 0.5 corresponds to an impulse one half orbit prior to impact). 
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Impact True Anomaly = 30 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 




Figure 18. Earth to Asteroid Separation 

2. Relative Maxima 

Impulse times ranging from 0 to 114 asteroid orbits prior to impact yield two 
maxima points as shown in Figure 1 8. Inspection of these maxima points reveal that they 
occur at the time of perihelion passage for the asteroid approximately one orbit prior to 
impact for impulse directions parallel and anti-parallel to the orbital velocity vector. This 
relationship holds true for all cases considered. 

3. Non-perihelion Impulse Direction 

Closer inspection of Figure 18 reveals that for impulse times less than that 
corresponding to one orbit prior to impact, the maximum separation is achieved if the 
impulse direction is not aligned either parallel or anti-parallel to the velocity vector. This 
effect is better shown in Figure 19, which is a contour plot of Figure 1 8. 



38 



1.5 



Impact True Anomaly - 30 deg, Orbital Eccentricity = 0.5, dv = 1 m/s 




A distinct shift in the direction of impulse for maximum separation as the time of 
impulse becomes closer to impact can be seen. This shift in direction arises from the two 
ways in which to achieve the desired separatioa If sufficient time prior time prior to 
impact exists, changing the speed of the asteroid on its orbit will shift the phase of the 
asteroid with respect to the Earth and avoid the impact. If the time prior to impact is 
short, changing the direction of the asteroid laterally with respect to its approach to the 
Earth is necessary to avoid the impact. For times prior to impact between these the 
former and the latter, a tradeoff between changing the orbital speed and displacing the 
asteroid laterally on its orbital path exists and is the cause of the shift in the optimal 
impulse direction. 

4. Periodic Growth 

Thus far, only impulse times within approximately one asteroid orbital period of 
impact have been considered. Extending the analysis to beyond this time period 
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demonstrates the manner in which the separation grows as a function of impulse time. 
Figure 20 shows this behavior for the same impact scenario as discussed above. Of note is 
the periodic growth over time and the peak displacements occurring at each perihelion 
point. This particular figure represents a slice of Figure 18 along the 0° impulse direction 
extended from 0 to 10 orbits prior to impact. 



Impact True Anomaly = 30 deg, Orbital Eccentricity = 0.5, dv = 1 m/s 




Figure 20. Periodic Growth of Separation Distance 

B. RESULTS 

The collection of the numerous impact scenarios modeled by the above method 
may be found in Appendix D. A study of these numerous model results yield the following 
general conclusion. 

1. Optimum Impulse Condition 

Assuming the optimum impulse condition is achievable in terms of a “real” mission 
sense (that is, the booster technology and energy delivery mechanism exist for asteroid 
mitigation), the optimum impulse point is located at the perihelion of the original asteroid 
orbit at least !4 orbit prior to impact. If time prior to impact permits, impulse at perihelia 
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multiple orbits prior have an even greater effect on separating the asteroid from the Earth 
at the given impact time. However, impact prediction becomes a problem in such cases as 
the validity of orbital prediction models in a general n-body problem comes into question. 

2. Optimum Impulse Direction 

If, due to the late time of detection, it is not possible to achieve even the first 
relative maximum for the optimum impulse condition, then there exists an optimum 
impulse direction at the time of impulse that maximizes the Earth to object separation. For 
fractions of an orbit from 0.2 < (t/P) < 0.9 there is a shift in the direction of impulse 
toward the orbital inward normal (for impulses that increase the orbital speed) and toward 
the orbital outward normal (for impulses that decrease the orbital speed) that yields a 
maximum in separation achievable. 

3. Time of Arrival Consideration 

Given the periodic nature of the impact problem, there exist conditions where it 
may be beneficial to delay a deflecting impulse until an optimal impulse condition occurs. 
For the scenario corresponding to Figure 20, if time permits delivery of an impulse two 
and one half asteroid orbits prior to impact, it would prove more advantageous to delay 
the impulse until only two orbits prior to impact in order to maximize the effect of the 
impulse. 

4. Detection Consideration 

The difficulty in realizing the use of an optimal impulse condition is that the 
detection of a colliding object may occur too late to achieve the most desirable condition. 
The earlier the detection the better the chance of deflection with a much smaller imparted 
energy. Unfortunately, the current search programs and record of detection have been 
yielding very short response times for NEO's having very close approaches. The current 
range of times for detection has been on the order of five days prior to closest approach to 
detection only after Earth passage. 
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VIII. APPLICATION TO A REAL CASE 



A. IMPULSE ACHIEVABLE 

The analysis of the energy coupling between an explosive yield and an asteroid 
performed by Ahrens and Harris (1994) shows that it may be feasible to deflect a NEO 
with an impulse magnitude from 1 cm s' 1 to a few m s' 1 for a globally threatening object of 
about 1 km in diameter. This impulse may be achieved using one of several methods. 
However, due to the relatively short warning times considered in this thesis, a nuclear blast 
appears to be the most efficient energy delivery mechanism The blast may be a standoff 
detonation or surface detonation. 

For the standoff detonation, the required explosive yield, W, may be determined by 
the approximate expression: 

10 3 AvD 3 



from Ahrens and Harris (1994), which has been modified to express yield in kT of 
equivalent TNT. The impulse, Av, is expressed in m s' 1 and the asteroid diameter, D, is in 
km. The efficiency of neutron production, n, from the nuclear blast lies between 0.03 and 
0.3. The standoff blast efficiency factor, A, is taken for an optimum standoff distance of 
0.4 asteroid radii with an efficiency of approximately 0.3. An order of magnitude analysis 
of the above approximation is presented in Table 3. 



Impulse (ms* 1 ) 


0.1 km diameter 


1 km diameter 


10 km diameter 


0.01 


0.1-1 kT 


100-1000 kT 


100-1000 MT 


0.1 


1-10 kT 


1-10 MT 


1-10 GT 


1.0 


10-100 kT 


10-100 MT 


10-100 GT 


10.0 


100-1000 kT 


100-1000 MT 


100-1000 GT 



Table 3. Impulse and Diameter v. Standoff Explosive Yield 



For the surface detonation, the required explosive yield, W, may be determined 
from the approximate expression: 

W =4xl0“ 9 AvM NEO , 
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also from Ahrens and Harris (1994) which again has been modified to express yield in kT 
of equivalent TNT. The impulse, Av, is expressed in m s' 1 and the asteroid mass, Mneo, is 
in kg. Aji order of magnitude analysis of the above approximation is presented in Table 4. 



Impulse (m s' 1 ) 


0.1 km diameter 


1 km diameter 


1 0 km diameter 


0.01 


60 T 


60 kT 


60 MT 


0.1 


600 T 


600 kT 


600 MT 


1.0 


6 kT 


6 MT 


6 GT 


10.0 


60 kT 


60 MT 


60 GT 



Table 4. Impulse and Diameter v. Surface Explosive Yield 



B. TOUTATIS 

Asteroid 4179 Toutatis will have multiple close approaches with the Earth over the 
next 15 years. It is of interest to apply the above analysis and methodology to Toutatis as 
if it were going to impact the Earth. 

To perform this analysis it must be assumed that the orbit of Toutatis lies in the 
ecliptic plane. This is not far from the true geometry of Toutatis' orbit where the orbital 
inclination is 0.47° out of the ecliptic plane. From the catalog of NEA orbital elements 
listed compiled by Tholen (1995), the semi-major axis and eccentricity of Toutatis are 
currently listed as 2.5154 and 0.6361, respectively. Assuming that Toutatis and the Earth 
will collide and that the Earth is in a circular orbit at 1 AU yields an impact at ±38.53° 
with respect to the perihelion of Toutatis. This is determined by solving Equation (3) for 
true anomaly. Figure 21 depicts this relative orientation of Toutatis with respect to the 
Earth. For the following analysis, the +38.53° impact location is chosen for modeling 
purposes. 
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Impact True Anomaly = 38.53 deg, Orbital Eccentricity = 0.6361 




Modeling the Toutatis collision in the manner described above yields the results 
displayed in Figure 22. From the JPL and NASA Photo Caption (1993) and Press Release 
(1996), the size of Toutatis is estimated to be approximately that of two attached spheres 
having diameters of about 4 km and 2.5 km. If both masses are combined, the effective 
spherical diameter is approximately 4.3 km. Using a mean asteroid density of 3000 kg m' 3 
results in a mass for Toutatis near 1.25xl0 14 kg. From the previous impulse analysis, an 
explosive yield of about 5 MT is required in the case of a surface detonation and an 
explosive yield from 9 to 90 MT, depending on neutron production, is needed for a 
standoff detonation to achieve a 1 cm s' 1 change in orbital speed. 

Using the 1 cm s' 1 orbital speed change determined above, a maximum separation 
of 1.64 Earth radii is the result of an impulse delivered 1.02 orbits prior to impact 
(perihelion passage of the prior orbit). While the separation is sufficient to cause the 
Toutatis to miss the Earth in this scenario, a greater margin of safety would be desirable. 

A larger separation may be achieved by either delivering an impulse greater than 1 cm s' 1 
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or delivering the impulse at an earlier perihelion passage. Recognizing the extremely large 
explosive yield requirements for increased impulse magnitudes, it would appear preferable 
to deliver the inpulse at an earlier time. This type of analysis demonstrates the necessity 
for detection of threatening asteroids many orbits prior to impact. 



Impact True Anomaly = 38.53 deg 
Orbital Eccentricity = 0.6361 
dv= 0 01 m/s 




Figure 22. Impulse Effects on Toutatis-Earth Impact 
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IX. FURTHER CONSIDERATIONS 



The above analysis shows promise as a tool for rapidly evaluating numerous 
scenarios for deflecting an asteroid that is going to impact the Earth. The possibility for 
further investigation utilizing this method presents itself in analyzing longer response times 
prior to impact, generalizing the method to a three-dimensional method and eccentricities 
other than l A. Additionally, more work is needed in analysis of very short response time 
impulse effects. 

A. LONG RESPONSE TIME 

The analysis presented above and in Appendix D have been performed primarily 
for impulse times between 0 and IV 2 orbits prior to impact. A few models were made for 
one impulse direction at times ranging from 0 to 10 orbits prior to impact. However, this 
investigation needs to be pursued further in search of general trends other than maximum 
separations occurring at perihelion points. 

B. THREE-DIMENSIONAL ANALYSIS 

The current model and method apply only to two-dimensional scenarios. It is of 
interest and merit to further generalize the analysis to the three-dimensional case. This will 
allow for orbital inclinations out of the ecliptic plane and better simulate a variety of real 
scenarios. 

C. ECCENTRICITY VARIATIONS 

In the above analysis, other than for Toutatis, the orbital eccentricities have been 
maintained at V 2 . An investigation into the effects of more circular orbits and highly 
eccentric orbits is in order. 

D. SHORT RESPONSE TIME 

The method presented is derived from a two-body representation of a more 
complicated physical system. This method is not valid for very short response times when 
the impacting object is within the Earth's sphere of influence. A further investigation is 
desirable to assess the effects of impulses in the three-body problem that arises when the 
object is detected close to the Earth. 
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APPENDIX A. DETAILED SOLUTION METHOD 



Given an Earth-asteroid collision with the following properties: 

A circular Earth orbit with semi-major axis (radius), a £ = 1AU 
(orbital eccentricity, e £ =0); 
an asteroid orbit with orbital eccentricity e v ; 
the true anomaly at time of impact, v topact U = 



and the time of the perturbing impulse, Av , prior to impact 



o 



V p i 






Av 



Find the minimum separation of the asteroid from the Earth in the vicinity of the old impact point. 
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Nomenclature 



Variables: 

a 

E 

e 

h 

F 

N 

n 

P 

P 

R 

r 

v 

x 

y 

Av 

a 

P 

v 



P 

Av 




semi-major axis 
eccentric anomaly 
orbital eccentricity 
specific angular momentum 
fraction 
number 

orbital mean motion 
orbital period 
orbital parameter 

asteroid to Earth separation distance, km 

orbit to Sun radial distance 

velocity 

x-position 

y-position 

impulse quantity 

angle w.r.t. x-coordinate axis 

gravitational mass parameter 

true anomaly 

asteroid to Earth separation distance, Earth radii 
shift in periapsis direction 

fraction of orbit from periapsis 
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Subscripts: 



days 

E 



inc 

impact 

min 

orbit 

P 

P 

pos 

Range 

sep 

Sun 

U 

x 

y 

Av 

ii 

± 

Symbols: 

A 

5 



number of days 

Earth 

i* element 

increment 

impact position 

minimum 

orbit count 

perturbed asteroid orbit 

periapsis 

position 

set of all i’s 

separation 

Sun parameter 

unperturbed asteroid orbit 

x-component 

y-component 

impulse 

parallel component 
normal component 

incremental amount 
incremental amount 
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Elliptical Orbit Relations 



Mean motion: 



n = 



Period: 



P = 



2k 



n 



Radial distance to gravitational focus (general): 

a(l -e 2 ) 

r = — — 

1 + ecosv 

Periapsis radius 

r a (l-f ecos v a ) 



r r = 



1 + e 



Semi-major axis: 



a = 

1 - e 

Eccentric anomaly as a function of true anomaly: 

e + cos v 



cosE = 



1 + ecos v 

True anomaly as a function of eccentric anomaly: 

cos E - e 

cosv = - ^ 

1 -ecosE 

Time since periapsis: 

E -esinE 



t = 



True anomaly: 



n 
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re e 

Position at a point in perifocal coordinates: 
? = (r cos v)x + (r sin v)y 
Velocity at a point in perifocal coordinates: 

v = 

Parameter of an orbit: 




Specific angular momentum: 

h=f xv 

Eccentricity vector in perifocal coordinates: 

1 

e = — 



|v| 2 — j^Jr-(r-v)v 



sin v)x + (e + cos v)y] 



Constants 

1AU = 1.4959787xl0 8 km, astronomical unit 
G = 6.67259x10 ” 20 km 3 kg -1 S “ 2 , gravitational constant 
M Sun = 1.9891x10' 30 kg , mass of the Sun 

ft Sun = GM Sun = 1.327124399355xl0 u km 3 s -2 , Sun gravitational parameter 
R e = 6.37814xl0 3 km , Earth radius 
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Detailed Solution Method 



I. Determine unperturbed asteroid orbital elements, 

a. Perihelion radius: 



r P ,u ~ 



a E (l + e u cosv impactu ) 



1 + e 



, km 



u 



where r tapact #u = a E 



b. Semi-major axis: 



PAJ . 

^TT “ ^ y km 



u 1-e, 



c. Mean motion: 



n u =■ 



I M' Sun 



, rad s’ ] 



d. Period: 



P TT =■ 



2k 



n T 



%s 



II. Determine the conditions at impact with respect to the unperturbed asteroid orbit, 
a. Eccentric anomaly at impact: 



( 



^ impact,U 



-1 



e U + COS V impact^ 



+ e u ^OS ^impact,U J 



(cos' 1 principle values are 0 < 0 < n ) 
b. Fraction of orbit since (or prior to) perihelion of impact: 









impact, p 



^ impact ,U ^ impact,U 

2n 



^ V impact, U < ^ > | r> 



<0 



u / • 

u impact ,p 
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III. Determine the conditions at impulse with respect to the unperturbed orbit. 



a. Fraction of orbit prior to (or since) perihelion of impulse: 






vPuJ 






vl uA 






Av,p u impact ,p u Av,impact 

b. Eccentric anomaly at impulse (solve via Newton-Raphson iterative method): 






pj 



Av,p 



E Av,U ~ e U COSE Av,U 

271 



c. True anomaly at impulse measured with respect to unperturbed perifocal coordinates: 



V Av,U =cos 



-l 



c °sE A vU e u 



1 -e v cosE 



Av,U J 



(cos -1 principle values are 0 < 0 < 7t ) 



L u' 



= N + F 

orbit ^ 1 orbit 



Av,p 



if 



~ ^ < E orbit < » ^ en V Av,U = orbit — V Av,U 




1 

~ 2 ~ E orbit ^ ^ * ^ en V Av,U = orbit ~ V Av,U 




> f 0 ^ F orbit < 2 > then V Av,U = 2 ^ N orbit + V AV/U 
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1 

if 2 - F °rbit < 1 » then 



V Av,U — 2 ti(N orblt +l) V Av,U 




d. Distance to focus (Sun) from impulse location: 

a u(l -e u) 



Av,l 



, km 



v ' u " (■ i \ 

(1 + e u cos V Av U J 

e. Velocity components at impulse with respect to unperturbed perifocal coordinates: 



^Su 



x,Av,U 



l u(l -e u) 



sin V Av,U > k™ s ’’ 



""• u ■jai|^ej)^ u + COSVav - u ) ,kmS " 



f. Velocity direction with respect to unperturbed perifocal coordinates 



\ 



a = tart 



. y,Av,U 
\ ^x,Av,U / 



(achieved utilizing “atan2(y,x)” numerical routine) 
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IV. Now perturb the asteroid with respect to velocity vector direction. 



a. Impulse velocity components with respect to unperturbed perifocal coordinates: 

Av x = AVu cos a - Av x sin a , km s'* 

Av y = Av B sin a + Av ± cos a , km s ' 1 
where, Av^ is in the direction of the velocity vector 

and, Av x is normal to the velocity vector in a right hand sense 

b. Velocity components after perturbation with respect to unperturbed perifocal coordinates 





y 



Av ± cos a 



Av ± sin a 




Av„ cos a 



v 
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V. Determine perturbed asteroid orbital elements from f and v . 

a. Position vector with respect to unperturbed perifocal coordinates: 

r Av,P = r Av,U = r x,Av,U^ + r y,Av,U Y + 

where, r X/AV/U = r Av ,u cos v av,u 

r y,Av,U = r Av,U V Av,U 

b. Velocity vector with respect to unperturbed perifocal coordinates: 

V 4v,P = V *,Av,pX + V y,4v,py + Oz , km s' 1 

c. Specific angular momentum vector with respect to unperturbed perifocal coordinates 

h p = r Av p x v Av p , km 2 s 1 

d. Perturbed orbit parameter: 




where |h p | 2 = h p *h p 



e. Perturbed eccentricity vector with respect to unperturbed perifocal coordinates: 

V \ 



e p = 



M'Su 



M'Su 



Av,P' 



V 



I r Av,pU 



r Av,P ( r Av,P * V Av,P ) V Av,F 



where |v AvP | 2 = v AvP * v Av#p , 

I r Av,pl = r Av,U 



f. Perturbed orbital eccentricity: 



e p ” V e P * e P 



g. Perturbed orbit semi-major axis: 
Pp 



a p=—— j^lcm 



h. Perturbed orbit mean motion: 
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n p = 



/M'Sun , .i 

[ — — , rad s 



i. Perturbed orbital period: 

_ 2n 
P p = , s 



VI. Determine the angle between the perturbed and unperturbed orbital eccentricity vectors. 

a. Unperturbed orbital eccentricity vector: 

e v =e u x + 0y + 0z 

b. Rotation of eccentricity vector due to impulse: 



Av = cos 1 



^ e p- e u A 

6p6;j ) 
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VII. Determine the impulse conditions with respect to the perturbed orbit, 
a. True anomaly at impulse with respect to the perturbed orbit: 



v 



Av,P 



= COS ' 1 



\ F (l-e 2 p ) 




T P 

V Av,P c P 


e P> 



(cos" 1 principle values are 0 < 0 < 7t ) 

if-1 <F orWt <--,then v Av<P = 27 i(N orbit -l)+v AvP 




if - 2 - F orbit < 0 > then v AV/P = 2%N olhit - v Av , P 




if 0 < F orbit < - , then v AV/P = 27tN orbit + v AV/P 




if “ £ F orbit < 1 , then v Av p = 2rt(N orbit + 1) - v 



Av,P 
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b. Eccentric anomaly at impulse with respect to perturbed orbit: 



E 



Av,P 



= cos -1 



" e p +cosv Av P ^ 
Vl + e p cosv AvP J 



(cos' 1 principle values are 0 < 0 < n ) 

c. Fraction of orbit of prior to (or since) perihelion of impulse with respect to perturbed orbit: 



f— 

^P P ; 



Av,p 



E Av,P~ e P SinE Av,P 

2 n 



if-l <F orbit < — ~ , then 



vP ) (N orbit 1 ) + ( P 

r Av,p r 



Av,p 





if 0 < F, 



orbit < 2 ’ then 






vP P y 



= N orbit + 






Av,p 



\P P y 



Av,p 




if 2 " F orbit < 1 > then 




= K*. 





Av,p 
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VIII. Determine a range of times, and hence positions, of the unperturbed asteroid in the vicinity of the 
impact point on the unperturbed orbit for mapping to corresponding times, and positions, of the perturbed 
asteroid on the perturbed orbit. 

a. Set time period: 

N days = ±2, days 

b. Set number of positions to map during time period: 

n pos = 201 

c. Establish limits of time interval: 






V*u/ 



^ days ^ f 24hrsV 3600s 

, Pu JlldayJv lhr 



d. Number of increments in time period: 



^ rnc ^ pos ^ 



e. Incremental time step: 






vPu/ 



= A 






Uu/ 



^ n in J 



f. Now have positions from -N^ to +N days in steps of 2N days /n mc parts of a day represented 
as fractions of a complete orbit. 






vP y y 



= —A 



Range 



vivHig,- +A 






vPu^ 



, in steps of 8 zr - 



g. Center this range of positions on the fraction of orbit since (or prior to) perihelion of impact: 
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r t ) 


r t > 




r t^i 






+ 




,pj 

u Range,impact 


l p J 


Range 


l P iJ 



impact, p 

h. Eccentric anomaly at each time coordinate (solve via Newton-Raphson iterative method): 






v p uA 



Range, impact 



E i,u ~ e u cosEju 
2n 



i. Determine the true anomalies at each time coordinate with respect to the unperturbed perifocal 
coordinates. 



Vtfj = cos 



-l 



cosE^-e^ 

1 — Gy COsEj u 



(cos' 1 principle values are 0 < 0 < 71 ) 



if 






. . <--,then v w =-27l + v i/U 

\ 1 U ' • r> • . ^ 

^ i. Range, impact 






< 0 , then Vj u = -v. 0 



^ i, Range, impact 



if 



> ~Z > then Vj u = 2;t - v i; 



VP,J 

u i,Range, impact 

j. Determine the radial distance to the focus (Sun) for these time coordinates: 
l u(l -e u) 



r * ,u l + e.-cosy,,, ’ 



km 



-U v i,U 

k. Determine the distance components at each time coordinate with respect to the 
unperturbed perifocal coordinates: 

x i,u = r i,u V i,U s ^ Cin 



Yi,u = r i,u sinv i/U9 km 
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IX. We now need to determine the true anomaly and fraction of orbit on the perturbed orbit that 
corresponds to the true anomaly and fraction of orbit at the impact position on the unperturbed orbit. 
Determine the positions and orbit fractions on the perturbed orbit that correspond to the n^ 
positions and orbit fractions of the unperturbed orbit. 

a. Fraction of orbit since (or prior to) perihelion of the “impact” position with respect to the 
perturbed orbit: 



t 



P y 



aJ av U, 



P p .impact ' v ' Av,p 

b. Establish the n^ orbital positions and fractions of orbit on the perturbed orbit that correspond 
exactly to the n^ positions of the unperturbed orbit: 

f 

jv 






Up/ 



-1 (* 

Ip p . 






vP P / 



Range, impact ' u ' Range r ' r ' impact.p 

c. Eccentric anomaly at each time coordinate (solve via Newton-Raphson iterative method): 



A 



P p J. p 

r i ,Range, impact 



E iP e p cosE iP 
2n 



d. Determine the true anomalies at each time coordinate with respect to the perturbed perifocal 
coordinates. 



cosEj p -e p 

V; p = COS — 

'' ll-e D cosE. 



i,P > 



(cos' 1 principle values are 0 < 0 < TC ) 




P i. Range, impact 



< — — , then V; p = — 2tc + v i<p 



if 






vl p/ . „ 

v 1, Range,: 



< 0, then v iP = -v iP 



Range .impact 
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if 



\PpJ . 



Range, impact 




then v iP = 2n — v i P 



e. Shift true anomalies from perturbed coordinates to unperturbed coordinates: 

V i,P,U = v i,p + Av 

e. Determine the radial distance to the focus (Sun) for these time coordinates: 



r D = 



,(l-e|) 



,p 1 + e p cos Vj p 



, km 



f. Determine the distance components at each time coordinate with respect to the 
unperturbed perifocal coordinates: 

Xi p = Tj p cos Vj p u , km 



y i,P “ r i,P S ^ n V i / P,U » ^ 



X. Now we need the positions of the Earth at the same n^ positions surrounding the impact position, 
a. True anomaly of the Earth at each of the n^ positions: 



\ 



^Range,E ^impact,E 



,Pu> 






Range 



b. Determine the distance components of the Earth at each time coordinate with respect 

to the unperturbed perifocal coordinates: 

Xj E = COS v i E ? km 

y iE =a E smViE ? km 
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XI. Now must handle the special case of 




For this case, set: 



X i,P = X i,U 




yi,p = yi,u 




f 1 1 


r t^| 


while < A 




> 

<3 





XII. Define the separation distances of the perturbed and unperturbed asteroid orbital positions from the 
Earth orbital positions. 

a. Perturbed asteroid to Earth separation distance: 



R 



sep,P 




, km 



b. Unperturbed asteroid to Earth separation distance: 




c. Scale the separation distances by the radius of the Earth: 



Psep,P 



Psep,U 



R 



sep,P 



R, 



Rsep,U 



Earth radii 



, Earth radii 
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XIII. Finally, select the minimum value over the orbital position range as the minimum separation 
distance. 

a. Minimum perturbed asteroid to Earth separation: 

( Pse P' P )min = ^(Psep^)’ Earth rad " 

b. Minimum unperturbed asteroid to Earth separation (By the definition of the problem this is 
zero, but this provides a good check of the method.): 

(psep.u)^ = m * rl (psep / u) ’ ra< ^* 
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APPENDIX B. ANALYTIC SOLUTION METHOD, MATLAB MODEL 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/ o o/ 0 o/ 0 o /o 

% Arbitrary Elliptical Orbit Intersecting Circular Earth Orbit 
% 

% MATLAB Script File Name: nwarborb.m 
% 

% Author: LT Jeffrey T. Elder 
% Naval Postgraduate School 

% November 1996 

% 

%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%%%%%%%%%%%%%%%% 



%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o%%% 0 /o 0 /o%%%%%%%%%%%%%%%%% 

% Initialize computational workspace 
clear 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o%%% 0 /o% 0 /o 0 /o% 0 /o%% 0 /o 0 /o% 0 /o 0 /o 

% Define constants 

AU = 1 .4959787e8; % astronomical unit 

musun = 1 .3271 24399355el 1; % gravitational parameter for the Sun 

epsilon = 5e-6; % Newton-Raphson tolerance 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%%% 

% Gather input data from user 

nuimpactU = input('Enter True Anomaly for Impact (deg): ’); 

eU = input('Enter Asteroid eccentricity: '); 

tPUdelv = input('Enter (t/P) of Impulse Prior to Impact: '); 

delV = input(’Enter the Magnitude of the Impulse dV, (m/s): *); 

thet = input('Enter the Direction of the Impulse wrt V, ccw+ (deg): '); 

dwi = delV*cos(thet*pi/180); 

dvni = delV*sin(thet*pi/180); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% O /o O /o 0 /o%%% 

% Allow for model performance evaluation 

%flops(0) 

%tic 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o%%% 0 /o 0 /o 0 /o 0 /o°/o 
% Set Earth orbital elements 
aE = 1.0* AU; 
nE = sqrt(musun/aE A 3); 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o%% 0 /o 0 /o%°/o%%% 
% Convert impact true anomaly to radians 
nuimpactU = nuimpactU*pi/180; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o% 0 /o 0 /o% 0 /o%%% 



% Determine unperturbed asteroid orbital elements 
rpU = aE*(l+eU*cos(nuimpactU))/(l+eU); 
aU = rpU/(l-eU); 
nU = sqrt(musun/aU A 3); 

PU = 2*pi/nU; 



% Perihelion radius 
% Semi-major axis 
% Mean motion 
% Period 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%%% 0 /o 0 /o 0 /o 

% Determine the conditions of impact wrt unperturbed orbit 

EimpactU =... % ... is line continuation 

acos((eU+cos(nuimpactU)). . . 

/(l+eU*cos(nuimpactU))); % Impact eccentric anomaly 

tPUpimpact =... 

(EimpactU-eU*sin(EimpactU))/(2*pi); % t/P of impact wrt perihelion 

if nuimpactU < 0 
tPUpimpact = -tPUpimpact; 
end 



%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%% 



tPUpimpulse = tPUpimpact-tPUdelv; 
Norbit = fix(tPUpimpulse); 

Forbit = tPUpimpulse - Norbit; 



% t/P of impulse wrt perihelion 
% Number of whole orbits 
% Fraction of orbits 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Allow for numerical simulation validation 

tao = tPUpimpulse*PU; % Start time for simulation 

tend = (tPUpimpact+O. 125)*PU; % End time for simulation 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o%% 0 /o%%%%% 0 /o%%% 0 /o% 0 /o%% 



% Determine true anomaly of impulse 
EimpulseU = pi/4; 
tptemp = Forbit; 
if tPUpimpulse < 0 
tptemp = -Forbit; 
end 

true = 0; 
while true = 0 

Iprime = 1 - eU*cos(EimpulseU); 
f = EimpulseU - eU*sin(EimpulseU) - 2*pi*tptemp; 
Elast = EimpulseU; 

EimpulseU = EimpulseU - f7fprime; 
if abs(EimpulseU-Elast) < epsilon, true=l;, end 
end 

if tPUpimpulse < 0, EimpulseU = -EimpulseU;, end 
nuimpulseU=... 
acos((cos(EimpulseU)-eU)./... 

( 1 -eU*cos(EimpulseU))); 



if (-1 < Forbit) & (Forbit < -0.5) 
nuimpulseU =... 

2*pi*(Norbit-l)+nuimpulseU; 
elseif (-0.5 <= Forbit) & (Forbit < 0) 
nuimpulseU = 2*pi*Norbit - nuimpulseU; 
elseif (0 <= Forbit) & (Forbit < 0.5) 
nuimpulseU = 2*pi*Norbit + nuimpulseU; 
elseif (0.5 <= Forbit) & (Forbit < 1) 
nuimpulseU = 2*pi*(Norbit+l) - nuimpulseU; 
end 



% Guess eccentric anomaly 
% allow for +/- tPUpimpulse 



% Newton-Raphson iteration 



% End N-R iteration 



% True anomaly at impulse 
% If t/P of impulse exceeds 
% one orbit, determine true 
% anomaly for multiple orbits 
% Conditions to properly 
% locate true anomaly 
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% Determine orbital distances at time of impulse 
rimpulseU = aU*(l-eU A 2)/(l+eU*cos(nuimpulseU)); 
ximpulseU = rimpulseU*cos(nuimpulseU); 
yimpulseU = rimpulseU*sin(nuimpulseU); 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o/oO/oO/ 0 o/ 0 o/ 0 o/ 0 o /o o /o o /o o /o 

% Perturb the asteroid 



pU = aU*(l-eU*eU); 
vximpulseU =... 

-sqrt(musun/pU)*sin(nuimpulseU); 
vyimpulseU =... 

sqrt(musun/pU)*(eU+cos(nuimpulseU)); 
dw = dwi/ 1000.0; 
dvn = dvni/1000.0; 

alpha = atan2(vyimpulseU, vximpulseU); 
vximpulseP =... 

vximpulseU+dvv*cos(alpha)-dvn*sin(alpha); 
vyimpulseP =... 

vyimpulseU + dw*sin(alpha) + dvn*cos(alpha); 



% Parameter of orbit 
% Unperturbed velocities 

% Impulse parallel to V 
% Impulse normal to V 
% Angle of V wrt x-axis 

% Perturbed velocities 



%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%%%% 0 /o%%%%%% 
% Determine perturbed asteroid orbit elements from R and V 

rimpulsePvec = [ximpulseU yimpulseU 0]; % Position vector 

vimpulsePvec = [vximpulseP vyimpulseP 0]; % Velocity vector 

hpvec = cross(rimpulsePvec, vimpulsePvec); % Angular momentum vec. 

pP = hpvec*hpvec'/musun; % Parameter of orbit 

ePvec = ( (vimpulsePvec* vimpulsePvec'-musun/rimpulseU)*rimpulsePvec. . . 

- (rimpulsePvec* vimpulsePvec')* vimpulsePvec )/musun; % Eccentricity vector 

eP2 = ePvec*ePvec’; 

eP = sqrt(eP2); % Eccentricity 

aP = pP/(l-eP2); % Semi-major axis 

nP = sqrt(musun/aP A 3); % Mean motion 

PP = 2*pi/nP; % Period 



%%%%%%%%%%%%%%%%%%%%%V0%%%%%%%%%%%%%%%%%%%%% 0 /0°/0%%%%%% 
% Determine angle perturbed orbit perihelion makes wrt unperturbed orbit perihelion 
dnu = acos( (ePvec*[eU 0 0] , )/(eP*eU) ); 
if ePvec(2) < 0, dnu = -dnu;, end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o°/o 
% Determine impulse true anomaly wrt perturbed orbit 
cnuimpulseP = aP*(l-eP2)/(rimpulseU*eP) - 1/eP; 

nuimpulseP = acos(cnuimpulseP); % True anomaly at impulse 

% If t/P of impulse exceeds 
% one orbit, determine true 
% anomaly for multiple orbits 

if (-1 < Forbit) & (Forbit < -0.5) % Conditions to properly 

nuimpulseP =... 

2*pi*(Norbit-l)+nuimpulseP; % locate true anomaly 

elseif (-0.5 <= Forbit) & (Forbit < 0) 
nuimpulseP = 2*pi*Norbit - nuimpulseP; 
elseif (0 <= Forbit) & (Forbit < 0.5) 
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nuimpulseP = 2*pi*Norbit + nuimpulseP; 
elseif (0.5 <= Forbit) & (Forbit < 1) 
nuimpulseP = 2*pi*(Norbit+l) - nuimpulseP; 
end 



%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o% 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 



% Determine impulse time wrt perturbed orbit 
EimpulseP = acos( (eP+cnuimpulseP)/(l+eP*cnuimpulseP) ); 
tPPpimpulse =... 

(EimpulseP-eP*sin(EimpulseP))/2/pi; 



if (-1 < Forbit) & (Forbit < -0.5) 
tPPpimpulse = (Norbit- l)+tPPpimpulse; 
elseif (-0.5 <= Forbit) & (Forbit < 0) 
tPPpimpulse = Norbit - tPPpimpulse; 
elseif (0 <= Forbit) & (Forbit < 0.5) 
tPPpimpulse = Norbit + tPPpimpulse; 
elseif (0.5 <= Forbit) & (Forbit < 1) 
tPPpimpulse = (Norbit+1) - tPPpimpulse; 
end 



% t/P of impulse 
% If t/P of impulse exceeds 
% one orbit, determine t/P of 
% impulse for multiple orbits 
% Conditions to properly 
% determine t/P of impulse 



%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o%% 0 /o 0 /o 0 /o 0 /o%% 0 /o%% 0 /o 0 /o 



ndaycoef = 1.5; 

Ndays = ndaycoef*tPUdelv*delV; 



loop = 2; 
while(loop >= 1) 

npos = max([201 fix(100*Ndays+l)]); 



% npos & Ndays define initial 
% mesh size 

% Make two passes to ensure 
% 'mesh' is properly sized 
% to give accurate minimum 



nuU=[]; 

nuP=[]; 



% Initialize true anomaly 
% matrices 



%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o% 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o%%%%%% 0 /o 0 /o 



% Create a mesh of orbit positions about the impact point on the 
% unperturbed orbit 
deltPU = Ndays*24*3600/PU; 
tPURangeimpact =... 

linspace(-deltPU,deltPU,npos)+tPUpimpact;% located in t/P space 
for i = 1 :npos 
true = 0; 

EU = tPURangeimpact(i)*pi; 
tptemp = tPURangeimpact(i); 

if tPURangeimpact(i) < 0, tptemp = -tPURangeimpact(i);, end 
while true = 0 
fprime = 1 - eU*cos(EU); 
f = EU - eU*sin(EU) - 2*pi*tptemp; 

Elast = EU; 

EU = EU - f/fprime; 

if abs(EU-Elast) < epsilon, true = 1;, end 
end 

if tPURangeimpact(i) < 0, EU = -EU;, end 
nuU(i) =... 



% Unperturbed mesh half width 
% Unperturbed mesh properly 

% Determine true anomaly of 
% unperturbed mesh points 



% Newton-Raphson iteration 



% End N-R iteration 
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acx>s((cos(EU)-eU)./( 1 -eU*cos(EU))); % True anomalies of 

% unperturbed mesh points 

if tP URangeimpact(i) < -0.5 % Ensure true anomalies are 

nuU(i) = -2*pi + nuU(i); % properly located 

elseif tPURangeimpact(i) < 0 
nuU(i) = -nuU(i); 
elseif tPURangeimpact(i) > 0.5 
nuU(i) = 2*pi - nuU(i); 
end 
end 

rU = aU*(l-eU A 2) ./ ( 1 +eU*cos(nuU)); % Positions of unperturbed 

xU - rU.*cos(nuU); % mesh points wrt unperturbed 

yU = rU.*sin(nuU); % coordinate frame 



%%%%%%%%%%%%%%%%%%%%%%%%%%%o/o%%0/ o o/ o o/ o o/ o o /o 0 /o o /o o/ o o /o o /o o /o o/ o o /o o /o o / 0 o /o o/ o o/ o o /o o /o 0 /o 



% Map the mesh of orbit positions about the impact point onto the perturbed orbit 


tPPpimpact = tPPpimpulse + tPUdelv*PU/PP; 


% t/P of impact wrt perihelion 


deltPP = Ndays*24 *3 600/PP; 


% Perturbed mesh half width 


tPPRangeimpact =... 


% Perturbed mesh properly 


linspace(-deltPP,deltPP,npos)+tPPpimpact; 0 /o located in t/P space 
for i = 1 :npos 


% Determine true anomaly 


true = 0; 


% of perturbed mesh points 


EP = tPPRangeimpact(i)*pi; 
tptemp = tPPRangeimpact(i); 

if tPPRangeimpact(i) < 0, tptemp = -tPPRangeimpact(i);, end 




while true = 0 


% Newton-Raphson iteration 


fprime = 1 - eP*cos(EP); 
f = EP - eP*sin(EP) - 2*pi*tptemp; 




Elast = EP; 

EP = EP - frfprime; 

if abs(EP-Elast) < epsilon, true = 1 end 
end 


% End N-R iteration 


if tPPRangeimpact(i) < 0, EP = -EP;, end 
nuP(i) =... 

acos((cos(EP)-eP)./( 1 -eP*cos(EP))); 


% True anomalies of perturbed 


if tPPRangeimpact(i) < -0.5 


% mesh points 
% Ensure true anomalies are 


nuP(i) = -2*pi + nuP(i); 


% properly located 


elseif tPPRangeimpact(i) < 0 
nuP(i) = -nuP(i); 
elseif tPPRangeimpact(i) > 0.5 
nuP(i) = 2*pi - nuP(i); 
end 
end 

nuPU = nuP + dnu; 


% Perturbed mesh true 


rP = aP*(l-eP2) ./ (l+eP*cos(nuP»; 


% anomalies wrt unperturbed 
% coordinate frame 
% Positions of perturbed mesh 


xP = rP.*cos(nuPU); 


% points wrt unperturbed 


yP = rP.*sin(nuPU); 


% coordinate frame 



%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o%% 0 /o 0 /o 0 /o 0 /o 0 /o 

% Map the mesh of orbit positions about the impact point on the Earth’s orbit 

delnuE = Ndays*24*3600*nE; % Earth mesh half width 
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nuRangeE =... 

linspace(-delnuE,delnuE,npos)+nuimpactU; 

xE = AU*cos(nuRangeE); 
yE = AU*sin(nuRangeE); 



% Earth mesh properly located 
% in true anomaly space 
% Positions of Earth mesh 
% points wrt unperturbed 
% coordinate frame 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%%%%% 0 /o 

% handle special case of tPUdelv < deltPU 
if tPUdelv < deltPU 

[v,loc] = min(abs(tPPRangeimpact-tPUpimpact+tPUdelv)); 
xP(l :loc) = xU(l :loc); 
yP(l :loc) = yU(l :loc); 
end 



%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o% 



% Determine the Earth to asteroid separation in Earth radii 
rhosepU =... 

sqrt((xU-xE). A 2 + (yU-yE). A 2)/6378.14; 
rhosepP =... 

sqrt((xP-xE). A 2 + (yP-yE) A 2)/6378.14; 

[rhosepUmin,inu] = min(rhosepU); 

[rhosepPmin,inp] = min(rhosepP); 



% Unperturbed orbit 

% Perturbed orbit 
% Minimum unperturbed 
% separation (must be zero) 

% Minimum perturbed separation 



%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o%% 0 /o% 0 /o 0 /o 0 /o 



if Ndays “ ndaycoef*tPUdelv*delV 
Ndays= 1.25*abs(tPURangeimpact(inp)*PU/3600/24 ... 
-tPURangeimpact(inu) *PU/3 600/24) ; 
end 

loop = loop - 1; 
end 



% Refine mesh 



% Allow for loop exit 
% End while(test) loop 



%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o%% 0 /o% 

rhosepPmin % Display minimum separation 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o 0 /o% 0 /o%% 

%toc % Display model performance 

%flops 
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APPENDIX C. NUMERICAL VALIDATION MODEL 



SIMULINK MODEL 

M time | 
time 



Q- 

Clock 




Unperturbed Orbit 



^MATLAB, 

1 Function 
orb eom2 




Integrator 
Initial Conditions 
Vx + dvx 
Vy + dvy 
X 
Y 



Demux2 



Perturbed Orbit 



In order to run the SIMULINK numerical integration model, the analytic method 
model must first be run. The initial conditions for the integrators are specified by the 
working variables in the analytic model. 
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MATLAB SCRIPT SUPPORTING SIMULINK MODEL 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

% orbeom2d.m 
% 

% Equations of motion for 2d orbit 
% 

% Author: LT Jeffrey T. Elder 
% Naval Postgraduate School 

% November 1996 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 



function xdot=orb_eom(x) 

mu = 1.327124399355ell; 
AU= 1.49597870e8; 
r2 = x(3) A 2 + x(4) A 2; 
mur2 = mu/r2; 
r = sqrt(r2); 

xdot(l) = -mur2*x(3)/r; 
xdot(2) = -mur2*x(4)/r; 
xdot(3) = x(l); 
xdot(4) = x(2); 



% Mass parameter for Sun 
% Astronomical unit 
% Radial distance squared 

% Radial distance 

% vx 1st ODE 
% vy 1st ODE 
%x 1st ODE 
%y 1st ODE 
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APPENDIX D. ANALYTIC METHOD SCENARIO MODELS 



The following collection of figures represents the analysis performed using the 
numerical routine developed to perform the analytic solution. The scenarios are evaluated 
for impact true anomalies from -180° to +180° referenced to the asteroid orbit perihelion. 
The "steps" between impact scenarios is generally in 30° increments, however, for clarity 
purposes, the step size is reduced to 10° or even 1° intervals at times. The primary 30° 
increment figures occur where a "surface plot" and impact scenario plot appear together. 
At other increment values, two surface plots appear together. 

For the surface plots at impact scenarios approaching ±180°, the surface appears 
rough due to the coarseness of the step size in impulse time and impulse direction. At 
smaller step sizes the behavior is smooth and well behaved. 
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Separation (Earth Radii) 



Impact True Anomaly = -180 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 




Impact True Anomaly = -180 Orbital Eccentricity = 0.5 
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Y (AU) Separation (Earth Radii) 



Impact True Anomaly = -150 deg 
Orbital Eccentricity * 0.5 
dv = 1 m/s 




Impact True Anomaly = -150 Orbital Eccentricity = 0.5 




X (AU) 
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Impact True Anomaly = -120 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 




360 



0*0^ C 



Impact True Anomaly = -120 Orbital Eccentricity = 0.5 




-1.5 -1 -0.5 0 0.5 1 

X(AU) 
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y (/\U) Separation (Earth Radii) 



Impact True Anomaly = -90 deg 




360 



Impact True Anomaly = -90 Orbital Eccentricity = 0.5 

i i i i | j r 

Semi-major aws = 1.33 AU 




i 1 I I 1 I l 

-2.5 -2 -1.5 -1 -0-5 0 0.5 1 

X (AU) 
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<s 

or 



o 

33 

2 



Impact True Anomaly = -60 deg 




360 



Impact True Anomaly = -60 Orbital Eccentricity = 0.5 




-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 

X(AU) 



88 



Y (AU) Separation (Earth Radii) 



impact True Anomaly = -30 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 




360 






Impact True Anomaly = -30 Orbital Eccentricity = 0.5 
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Separation (Earth Radii) Separation (Earth Radii) 



Impact True Anomaly = -10 deg 
Orbital Eccentricity = 05 
dv = 1 m/s 




Impact True Anomaly = -5 deg 
Orbital Eccentricty = 0.5 




360 
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Separation (Earth Rad.i) Separation (Earth Radii) 



Impact True Anomaly = -2 deg 




Impact True Anomaly = -1 deg 
Orbital Eccentricity = 0.5 
dv= 1 m/s 
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Y (AU) Separation (Earth Radii) 



impact T rue Anomaly = 0 deg 
Orbital Eccentricity = 0.5 
dv= 1 m/s 




Impact True Anomaly = 0 Orbital Eccentricity = 0.5 




X (All) 
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Separation (Earth Radii) Separation (Earth Radii) 




Impact T rue Anomaly = 2 deg 
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Impact T rue Anomaly = 5 deg 
Orbital Eccerrtntity = 0.5 
dv= 1 m/s 



.. n 




360 



Direct C 



Impact T rue Anomaly - 10 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 
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Impact True Anomaly = 30 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 




Impact True Anomaly = 30 Orbital Eccentricity = 0.5 
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y (All) Separation (Earth Radii) 



Impact True Anomaly = 60 deg 
Orbital Eccentricity = 0.5 
dv= 1 m/s 




Impact True Anomaly = 60 Orbital Eccentricity = 0.5 

I | | | 

Semi-major aws = 1.67 AU 



T 



nr 



t r 



1.5 



0.5- 



0 - 



-0.5- 



-1 



-1.5 ■ 




_i i_ 



-3 -2.5 



-1.5 



-1 -0.5 

X (AU) 



0.5 
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Y (AU) Separation (Earth Radii) 



Impact True Anomaly = 90 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 




1.5- 



i 1 1 

Semi-major axis = 1 .33 AU 



Impact True Anomaly = 90 Orbital Eccentricity = 0.5 

i i i i i i 



1 - 



0.5 



-0.5 



-1 




-1.5- 



_l L. 



i 



J L_ 



-2.5 -2 -1.5 



-0.5 
X (AU) 



0.5 1 
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Y (AU) Separation (Earth Radii) 



impact True Anomaly =120 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 




360 



Off **'*' 1 



Impact True Anomaly = 120 Orbital Eccentricity = 0.5 




-1.5 -1 -0.5 0 0.5 1 

X (AU) 
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Y (AU) Separation (Earth Radii) 



Impact T rue Anomaly = 1 50 deg 
Orbital Eccentricity = 0.5 




360 



Impact True Anomaly =150 Orbital Eccentricity = 0.5 
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Impact T rue Anomaly = 1 70 deg 
Orbital Eccentricity = 0.5 




360 



Impact True Anomaly = 173 deg 




360 
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Separation (Earth Radii) Separation (Earth Radii) 



Impact True Anomaly = 175 deg 
Ortwtal Eccentricity = 0.5 




Impact True Anomaly =178 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 
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Impact T rue Anomaly = 180 deg 
Orbital Eccentricity = 0.5 
dv = 1 m/s 




Impact True Anomaly =180 Orbital Eccentricity = 0.5 
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Impact True Anomaly = 0 deg, Orbital Eccentricity = 0.5, dv = 1 m/s 
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